Digits := 80:
with(Groebner):
read("models.txt"):
filtered := [op(models)]:
nops(models);
<Text-field style="Heading 1" layout="Heading 1">Code to find numerical solutions of bivariate zero dimensional system</Text-field>eqf := proc(a,b)
if not (type(a,numeric) and type(b,numeric)) then return false: fi:
return evalb(abs(a-b)<0.0001):
end:
inf := proc(a,S) return `or`(op(map(eqf,S,a))); end:
iszer := proc(a) if eqf(a,0) then return 0 else return a: fi: end:
zerpol := proc(p) local cfs,xs,t;
cfs := [coeffs(p,x,'t')]:
xs := [t]:
return add(iszer(cfs[k])*xs[k],k=1..nops(cfs));
end:
fsolvesys := proc(G) local pts,k,ys,j,tp,xs;
ys := {fsolve(op(remove(has,G,x)))};
pts := {}:
for k in ys do
xs := map(a->{fsolve(a)},remove(a->a=0,map(zerpol,subs(y=k,G))));
tp := xs[1];
for j in xs[2..-1] do
tp := select(inf,tp,j):
od:
pts := pts union {seq([j,k],j=tp)}:
od:
end:LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2I1EhRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYn
<Text-field style="Heading 1" layout="Heading 1">Find Min Polys of alpha</Text-field>minPs := proc(S)
local alpha,cycPol,chi,chi_x,chi_y,i,k,pos,GB0,GBc,ys,xs,c,p,check,mP,mC,G,Q,Qpol,P,sls,t,u;
# Define characteristic polynomial of S
chi := add(x^k[1]*y^k[2],k=S):
# Compute partials and their numerators
chi_x := numer(normal(diff(chi,x)));
chi_y := numer(normal(diff(chi,y)));
GB0 := remove(has,Basis([chi_x,chi_y,1-u*x*y],plex(u,x,y)),u):
sls := select(a->(a[1]>0) and (a[2]>0), fsolvesys(GB0)):
if nops(sls)>1 then error "Multiple positive solutions": fi:
if nops(sls)=0 then error "No positive solution!" fi:
pos := {x=sls[1][1],y=sls[1][2]}:
p := subs(op(pos),chi);
# Get minimal poly of p
G := Basis([chi_x,chi_y,numer(t-chi),u*x*y-1],lexdeg([u,x,y],[t])):
mP := factor(op(remove(has,G,{x,y,u}))):
# Find the min factor and check it annihilates p
check := 0:
if not irreduc(mP) then
check := -1:
for i in factor(mP) do
if abs(subs(t=p,i))<0.00001 then check := 0: mP := i: break: fi:
od:
fi:
if check=-1 then error "Can't find min poly of mP" fi:
Q := simplify( diff(chi,x,y) / sqrt(diff(chi,x,x)*diff(chi,y,y)) ):
Qpol := simplify( diff(chi,x,y)^2 / (diff(chi,x,x)*diff(chi,y,y)) ):
P := numer(t^2-Qpol):
c := evalf(subs(op(pos),Q)):
GBc:=Groebner[Basis]([chi_x,chi_y,P,1-u*x*y*denom(Qpol)],lexdeg([x,y,u],[t])):
mC:=factor(op(remove(has,GBc,{x,y,u})));
# Find the min factor and check it annihilates c
check := 0:
if not irreduc(mC) then
check := -1:
for i in factor(mC) do
if abs(subs(t=c,i))<0.00001 then check := 0: mC := i: break: fi:
od:
fi:
if check=-1 then error "Can't find min poly of c" fi:
# Find alpha
alpha := evalf(-1-Pi/arccos(-c));
cycPol := expand(x^degree(mC,t)*subs(t=(x^2+1)/(2*x),mC)):
cycPol := mul(k[1],k=sqrfree(cycPol)[2]);
if not irreduc(cycPol) then error "R(x) not irreducible": fi;
return [p,primpart(mP),c,primpart(mC),alpha,primpart(cycPol)];
end:# Takes a long time to run the code (1 - 2 hours)
# We read in the results instead
read "allPols.txt":
nops(allPols);
#for k from 1 to 12883 do
# PP[k] := minPs(models[k])[6];
#od:allPols[123];
minPs(models[123])[6];# All the polynomials R(x) found are irreducible
nops(remove(irreduc,allPols));# The maximum degree is 72
max(map(degree,allPols));k := 'k':
# It is sufficient to compare the R(x)'s to the first 400 cyclotomic polys
# The following is a lower bound on the degrees of cyclotomic polynomials
ceil(subs(k=349,k/(exp(gamma)*log(log(k))+3/log(log(k)))));
with(numtheory):
CPs := [0,seq(cyclotomic(k,x),k=1..500)]:ratIndex := {}:
for k from 1 to nops(allPols) do
if allPols[k] in CPs then ratIndex := ratIndex union {k}: fi:
od:ratExp := map(a->models[a],ratIndex):nops(ratExp);# lcm of denominators is 30 -- use this to get exact fractions
for k in ratExp do
printf("%a, %a\134n",k,round(30*minPs(k)[5])/30);
od:# Go backwards to verify that the calculation is correct (independently verify c)
for ST in ratExp do
alp := round(30*minPs(ST)[5])/30;
chi := add(x^j[1]*y^j[2],j=ST);
cc := -simplify(cos(-Pi/(alp+1)));
pts := simplify(op(select(a->evalf(subs(a,x)>0) and evalf(subs(a,y)>0), map(allvalues,[(solve([diff(chi,x),diff(chi,y)]))]))));
evalf(subs(pts,diff(chi,x,y) / sqrt(diff(chi,x,x)*diff(chi,y,y)))-cc);
printf("%a, ",%);
od:LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2I1EhRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYn